Optimized Jastrow-Slater wave functions for ground and excited states: 
Application to the lowest states of ethene 
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A quantum Monte Carlo method is presented for determining multi-determinantal Jastrow-Slater 
wave functions for which the energy is stationary with respect to the simultaneous optimization 
of orbitals and configuration interaction coefficients. The approach is within the framework of the 
so-called energy fluctuation potential method which minimizes the energy in an iterative fashion 
based on Monte Carlo sampling and a fitting of the local energy fluctuations. The optimization of 
the orbitals is combined with the optimization of the configuration interaction coefficients through 
the use of additional single excitations to a set of external orbitals. A new set of orbitals is then 
obtained from the natural orbitals of this enlarged configuration interaction expansion. For excited 
states, the approach is extended to treat the average of several states within the same irreducible 
representation of the pointgroup of the molecule. The relationship of our optimization method with 
the stochastic reconfiguration technique by Sorella et al. is examined. Finally, the performance of 
our approach is illustrated with the lowest states of ethene, in particular with the difficult case of 
the l 1 Bi u state. 
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I. INTRODUCTION 

Quantum Monte Carlo (QMC) methods have been 
successfully employed over the last decade to compute 
ground state electronic properties of large molecules and 
solids 1 -. Compared to other electronic structure ap- 
proaches, QMC has the advantage that it can be applied 
to sufficiently large systems and still provide an accurate 
description of both dynamical and static electronic cor- 
relation. Recently, QMC methods have also been used 
for the calculations of excited statesS* 3 *^, a rather novel 
direction where relatively little experience exists. 

In both variational Monte Carlo (VMC) and diffusion 
Monte Carlo (DMC), the many-body trial wave function 
determines the quality of the calculation and is com- 
monly chosen to be of the Jastrow-Slater type: a single- 
or multi-determinant wave function is multiplied by a 
correlation Jastrow factor to partly account for dynamic 
electronic correlation. While the Jastrow factor is gen- 
erally optimized within QMC using the variance mini- 
mization method^, the determinantal component of the 
wave function is computed with methods such as Hartree- 
Fock, multi-configuration self consistent field (MCSCF, 
CASSCF) or small scale configuration interaction (CI), 
and, in most cases, is left unchanged when the Jastrow 
factor is added. However, the Jastrow factor is a positive 
function of the inter-particle coordinates and, therefore, 
does not change the nodal surface of the trial wave func- 
tion. Consequently, the fixed-node DMC energy is solely 
determined by the determinantal part of the trial wave 
function^. Any approach aimed at improving the fixed- 
node DMC energy using a Jastrow-Slater trial wave func- 
tion must therefore reoptimize the determinantal part in 
the presence of the Jastrow factor. 

Recently, the energy fluctuation potential (EFP) ap- 
proach to the optimization of Jastrow-Slater wave func- 
tions has been developed and applied to ground state cal- 
culations: the energy has been minimized with respect to 



a subset of the orbital parameters in single-determinant 
wave functions-, the CI coefficients in multi-determinant 
wave functions 9 , and the Jastrow parameters in periodic 
systems 1 -. In the EFP method, the optimization is per- 
formed iteratively through Monte Carlo sampling and 
the fitting of the energy fluctuations. Here, we extend 
the EFP method to the simultaneous full optimization 
of CI coefficients and single-particle orbitals by reformu- 
lating the optimization of the orbitals as a super-CI ap- 
proach^ for Jastrow-Slater wave function: the orbital 
variations are expressed as single excitations to a set of 
external orbitals and the improved orbitals are then ob- 
tained from the natural orbitals of this super-CI expan- 
sion. Moreover, for excited states, we further develop 
the EFP technique to treat an average of several states 
within the same irreducible representation of the point- 
group of the molecule. With these additions, the method 
closely resembles the state-averaged MCSCF approach^. 
When one treats multiple states, the resulting Jastrow- 
Slater wave functions are however not orthogonal since 
orthogonality is only ensured for the determinantal parts. 

To illustrate the performance of our optimiza- 
tion approach, we study the 1 Bi„ state of ethene 
which has been the subject of much theoretical 
debate and reconsideration over more than twenty 
yeara 13 i 14 i 15 i 16 i 17 i 18 i 19 i 20 i 21 i 22 i 23 . The difficulty in describ- 
ing the l 1 Bi„ can be mainly ascribed to the fact that any 
method insufficiently accounting for electron correlation 
will strongly mix this valence state with Rydberg states. 
Explicit inclusion of dynamic correlation in the reference 
wave function appears to be necessary to avoid such mix- 
ing but this renders the construction of a reference for 
highly correlated quantum chemical calculations partic- 
ularly complicated. In addition, the result is sensitive to 
the choice of the basis and to the optimization of the oc- 
cupied a- orbitals which cannot all be treated as frozen. 
QMC calculations for the l 1 Bi„ state are also affected by 
similar problems and, to obtain a satisfactory description 
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of this state, we find it necessary to optimize the deter- 
minantal part of the wave function in the presence of 
the Jastrow factor, that is to include the feedback of dy- 
namic correlation on the determinantal reference. The 
Rydberg character of the initial trial wave function is 
successfully corrected by our EFP optimization, and the 
DMC excitation energy is lowered by about 0.5 eV upon 
optimization of the orbitals and brought in very good 
agreement with the most sophisticated quantum chemical 
results. Our optimization approach is very efficient and 
robust, converging in a small number of iterations even 
when the initial trial wave function is worsened by mak- 
ing the single-particle basis more diffuse. We find that 
the EFP optimization has superior convergence proper- 
ties than the stochastic reconfiguration method^Si by 
Sorella et al. Finally, we consider two Rydberg-like states 
2 1 A 9 and 2 1 Bi„ as examples of states which are not the 
lowest ones in their irreducible representations and for 
which we successfully employ our extension of the EFP 
approach to the simultaneous optimization of multiple 
states. In a related paper—, we further investigate the 
state-average EFP approach by studying conformational 
changes in formaldimine, formaldehyde and a protonated 
Schiff-base model, and find that the optimization of ex- 
cited states in the absence of symmetry constraints al- 
ways works reliably. 

In Sec. [H] we review the EFP approach and introduce 
our method to fully optimize CI and orbital coefficients 
for ground and excited states. The relationship with the 
stochastic reconfiguration method by Sorella et al. is also 
examined. Computational details are given in Sec. IIIII 
and the numerical results for ethene are shown in Sec. II VI 
Finally, we present our conclusions in Sec. [Vj 



II. OPTIMIZATION OF JASTROW-SLATER 
WAVE FUNCTIONS 



The trial wave functions commonly used in quantum 
Monte Carlo calculations arc of the Jastrow-Slatcr form: 



*t = J $t = J 



i.Ci 



(1) 



where J is the Jastrow factor which explicitly depends on 
electron-electron separations and partially accounts for 
dynamic correlation. The CVs are a set of spin-adapted 
configuration state functions (CSF) expressed as linear 
combinations of Slater determinants: 



Ci 
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(2) 



For a given Jastrow factor, we want to minimize the en- 
ergy with respect to the CI coefficients Ci and the orbitals 
in the determinantal part of the wave function. 



A. Optimization of CI coefficients 

Let us assume that the determinantal component of the 
wave function <E>t (Eq. ^1 is an eigenstate of a CI Hamil- 
tonian in the basis of the configuration state functions 
C{. Since the configuration state functions Ci and the 
orthonormal eigenstates <I>i of this Hamiltonian span the 
same space, we can express the infinitesimal variations in 
the Ci coefficients of $t as variations with respect to the 
eigenstates $i other than <1>t: 



(3) 



fc>0 



where, without loss of generality, we have set <!>t equal to 
the lowest eigenstate $o- The corresponding derivatives 
of the energy are given by: 



dE 
d5 k 



<5=0 



8=0 



= ((El — E) (Ok — Ok)) 



(4) 



where (■) denotes the average with respect to the square 
of trial wave function, | \E r -x' | 2 , which can be conveniently 
computed by Monte Carlo sampling. We defined E — 
(E L ) and O k = (O k ) where 



O k — — and E L = — — . 
<j> \]/ T 



(5) 



The energy is stationary with respect to variations of the 
CI coefficients if these derivatives are zero. As shown in 
Eq. the fluctuations of the local energy become then 
uncorrelated with the fluctuations of the functions Ok, 
which means that the remaining fluctuations of the local 
energy cannot be further reduced by adding some com- 
bination of the functions Ok- 

The energy fluctuation potential (EFP) method is 
based on this last observation, and reformulates the en- 
ergy minimization problem as a least-squares fit of the 
fluctuations of the local energy with an arbitrary combi- 
nation of the functions Ok- 



X 2 = ((E L -Y^VkOk) 2 ) . 



(6) 



A minimization of \ 2 with respect to the parameters Vk 
leads to the following set of linear equations^: 



(E L O rn ) = J2 Vk ( 0k0 " 



(7) 



The equation for the trial state $o is simply 

E = V + J2 V k(Ok), (8) 



k>0 



and can be used to eliminate Vo from the other equations 
which, for m > 0, become: 



(AEAO m ) = V k (AO k AO m ) 



(9) 



fc>0 
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with AE = E L - {E L ) and AO m = O m - (O m ). The left- 
hand side of these equations correspond to the derivatives 
of the energy with respect to variations in the CI coef- 
ficients (Eq. 0}. Therefore, the fitting parameters Vfc>o 
are all zero if and only if all the derivatives of the energy 
are zero2&. 

For an arbitrary trial wave function \Pt = J&o, the 
parameters Vk which solve these linear set of equations 
will be different from zero. In order to understand how 
to use the coefficients Vk to obtain a new set of CI coeffi- 
cients (or Si coefficients in Eq. [3J| , let us first consider the 
case when there is no Jastrow factor and suppose that 
$o is the eigenstate <J > q°' > of an incorrect CI Hamiltonian 
W<°>: 

^=£4 0) i^°Wi, do) 

where the states $^°' ) are not the CI eigenstates for the 
given set of configuration state functions Ci . The coeffi- 
cients are easily obtained from Eq. ^\ as 

V^ = {^\H\^). (11) 

Therefore, the corrections V k ^ correspond to off-diagonal 
elements of the correct CI Hamiltonian coupling to the 
lowest state. With the use of these coefficients, a new 
Hamiltonian is constructed as 

H (n+i) =W (») + ^vW(|$W)($(«)| + / l , ( .) i (12) 

fc>0 

The updated Hamiltonian is diagonalized, yielding a new 
set of states and the procedure is repeated until 

convergence. An iterative approach is necessary because 
only the row and the column of the CI Hamiltonian cor- 
responding to $o are corrected at each step. At conver- 
gence, one obtains 14 = (^E>o l^l^fc) = and the correct 
CI eigenstate $o- Therefore, we have devised an itera- 
tive scheme to improve on the starting wave function and 
perform a CI calculation for $o- 

After reintroducing the Jastrow factor, we proceed mo- 
tivated by the above scheme developed in the absence 
of the Jastrow factor and by the observation that, if 
= J^o were an eigenstate of the Hamiltonian H, 
$o would be a right eigenstate of the operator J~ l TiJ 
with the same eigenvalue. Therefore, for Jastrow-Slater 
wave functions, the EFP method constructs iteratively an 
effective Hamiltonian 7i c ff which approximates J~ Xr H.J 
as far as the action on the trial state $o is concerned. 
One can then interpret Eq. as obtained from a least- 
squares fit of the local energy of the effective Hamilto- 
nian r T~LJ acting on the determinantal part of the 
trial wave function 

This interpretation allows us to use the corrections V k to 
update the approximate effective Hamiltonian TL c s in an 
iterative scheme as in the absence of the Jastrow factor. 



The initial guess for 7i c ff is arbitrary but, if we assume 
that an initial CI or MCSCF calculation yields a wave 
function not too far from the final one, we can construct 
a reasonable starting Hamiltonian using a complete set 
of states $fe (i.e. as many states as CSF) and the associ- 
ated energies E k from such a calculation. Typically, the 
QMC trial state <J>o is truncated according to a thresh- 
old on the CSF coefficients. Since the corrections to TC e s 
are sampled from J~$o, where $o is assumed to be an 
eigenstate of 7i c ff, it is important to keep H c s consistent 
with $o- To this end, 7i e ff is modified to yield the trun- 
cated $o as an eigenstate by setting to zero the matrix 
elements between the CSF's included in <!>o and the ones 
omitted. 

As the iterative process proceeds, the effective Hamil- 
tonian will contain more and more contributions sampled 
within quantum Monte Carlo and thus incorporate sta- 
tistical noise. Consequently, symmetries inherent to the 
true Hamiltonian could by broken. This can be avoided 
by using symmetry adapted expansion functions instead 
of simple Slater determinants. If molecular point groups 
are restricted to Abelian groups, individual determinants 
of symmetry adapted orbitals are proper spatial basis 
functions and only spin symmetry has to be considered. 

Recently, Umezawa and TsuneyukiSS proposed an opti- 
mization method for a single determinant Jastrow-Slater 
wave function based on the minimization of the quan- 
tity f dK 3N [n eS D - E gucss D) 2 where H cB = J^HJ. 
Within this approach, the derivatives of the energy with 
respect to the parameters in the determinant are in gen- 
eral not zero upon convergence. 

B. Optimization of Orbitals 

In order to obtain more compact wave functions, it is 
important to optimize the orbitals comprising the Slater 
determinants along with the coefficients in front of these 
determinants. Within the traditional quantum chem- 
istry framework, this approach is known as the MCSCF 
method, which is employed in the generation of reference 
wave functions for both complete-active-space second- 
order perturbation theory (CASPT2) and MRCI (multi- 
reference CI). Orbital optimization can be achieved ei- 
ther by using Lagrange multipliers, a method similar but 
much more involved then the SCF method, or using a 
so-called super-CI expansionii*2S. Within the latter ap- 
proach, the original reference space is augmented by all 
possible single excitations with respect to a set of exter- 
nal orbitals, leading to the super-CI wave function 




The notation C\ means that in the reference CSF C; 
an electron from the spin-orbital k is promoted to the 
external spin-orbital I. By virtue of the Brillouin-Levy- 
Berthier theorem^, it can be shown that the orbitals 
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of the reference space are optimal with respect to varia- 
tions in the external space if the coefficients in front of 
appropriate combinations of single excitations vanish. If 
convergence has not been reached, the coefficients cm can 
be used to generate improved reference orbitals. Here, we 
follow Rucdenberg et alm^ in generating new orbitals via 
natural orbitals of the super-CI wave function. 

We reformulate the minimization problem of the or- 
bitals in a Jastrow-Slater wave function as a super-CI 
approach. The resulting procedure consists of the follow- 
ing steps: 

Step 1. For a given trial wave function = J$a, ah 
single excitations out of the determinantal component $0 
are generated. In this enlarged space of CSF's, the super- 
CI Hamiltonian is set up and modified to yield $o (and 
not $sci) as an eigenstate. Thus, a starting effective 
super-CI Hamiltonian 7"4ci (Ea.llUp is obtained with $o 
as $q ^ . The quantities appearing in the linear equations 
(Eq. |5J) are then sampled in a variational Monte Carlo 
run with wave function \1/t, and the corrections Vk are 
obtained and added to the original Hamiltonian H S qj. 
The new effective Hamiltonian HgQj is diagonalized to 

yield a new state ^o^- 

Step 2. The natural orbitals (i.e. the eigenvectors of the 
single particle density matrix) are computed for the state 
The new reference orbitals are then obtained by 
demanding that the natural orbitals of the new reference 
wave function coincide with the natural orbitals of the 
super-CI wave function as explained in Appendix 1X1 

Step 3. After the orbital update is performed, the 
super-CI Hamiltonian needs to be recalculated in the ba- 
sis given by the new orbitals. While this is of course 
possible in the context of usual MCSCF, it is not in the 
case of the EFP since the effective Hamiltonian consid- 
ered here contains contributions sampled within QMC. 
Instead, we transform the effective Hamiltonian approxi- 
mately by projecting it onto a set of CSF consisting of the 
new orbitals but having the same occupation patterns as 
the old CSF. Subsequent QMC sampling of the matrix 
elements of the effective Hamiltonian as described above 
will correct this approximation. 

Let us denote with fa, Ci and Di the old orbitals, 
CSF's, and determinants, respectively, and with a tilde 
the corresponding new quantities. We project the old set 
of CSF's onto the new ones 



C, 



k 



Gik Cfc > 



(15) 



where the expansion coefficients are given by e^. = 
(Ci\Ck)- If the reference space is complete (CASSCF), 
this relation would be exact. Otherwise, it leads to a 
projection of the new effective Hamiltonian onto the new 
set of CSF's Ci. 

Since the CSF's are linear combinations of determi- 
nants (Eq. |2J) with the same set of coefficients dik for the 
old and new orbitals, the expansion coefficients eik can 



be evaluated as 



(16) 



The overlap of two Slater determinants is computed as 
the determinant of the overlap matrix of the spin orbitals 
occupied in the determinants under consideration: 



(Di\D k ) = det(0&) det(6f fe ) , 



(17) 



where Q? k = (4>f\<j)%), with fa* the i-th orbital of spin a 
in the determinant. 

Finally, the effective Hamiltonian is expressed in the 
projected basis Ci as 



(a|«eff|Q) = J2J2 eikejiiCklHenlQ] 



(18) 



To illustrate the effect of the projection, let us con- 
sider a system of two electrons and two orbitals fa , fa 
where the reference determinant ('active space') is fa^faf, 
and the single excitations into fa ('external space') yield 

With a new orbital 



the determinants faffa^ and 



'2 4>i- 



fa = c\fa + cifa-, the (projected) matrix element of the 
new reference determinant (f>i4>i is missing the contribu- 
tion from the determinant fa^fa^ which is proportional to 
c\ and is small if fa is close to fa. In general, if the new 
orbitals do not significantly differ from the old ones (i.e. 
the occupation numbers of the corresponding natural or- 
bitals are close to 0, 1, and 2, respectively), the projection 
approximation is expected to be quite accurate. 



C. Multiple states 

Even though, in the above description of the EFP 
method, we choose the trial state \Pt to be the lowest 
state 3<&q, this is not a requirement. Since the EFP 
conditions lead to a state for which the energy is station- 
ary with respect to parameter variations, arbitrary states 
could be optimized if a close enough starting guess is pro- 
vided. However, in optimizing a higher state, one some- 
times faces the so-called root flipping problem known 
from CI and MCSCF calculations of excited states: the 
higher state is lowered so much in energy that it ap- 
proaches (and mixes with) an initially lower state of the 
same symmetry which is not being optimized. When 
this happens, the procedure either converges to the lower 
state or, more likely, does not converge at all. 

Within the MCSCF framework, a well established solu- 
tion to the root flipping problem is the so-called state av- 
eraged MCSCF (SA-MCSCF) method 12 . In SA-MCSCF, 
the optimized quantity is the weighted average of the en- 
ergies of the states under consideration: 



E, 



(19) 
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where the weights Wi are fixed and Wj = 1. The 
multi-determinant wave functions ^ depend on their in- 
dividual sets of Cl-coefficients a k but on a common set 
of orbitals. Therefore, if the averaged energy E$a is sta- 
tionary with respect to all parameter variations, the in- 
dividual state energies E% are stationary with respect to 
variations of the Cl-coefficients but not with respect to 
variations of the orbitals. In SA-MCSCF, the wave func- 
tions are kept orthogonal and a generalized variational 
theorem applies^. 

The concept of state averaging is introduced into the 
EFP framework by considering several states € A, 
and associated weights Wi (^2 i Wi — 1) instead of a sin- 
gle trial state \&t- ln the EFP procedure, one needs to 
modify how the corrections to the effective Hamiltonian 
are obtained and how the natural orbitals are computed. 

For multiple states, we start from a least squares prob- 
lem like in Eq. but with the single trial function re- 
placed by the ensemble ^ — 3<&i,i G A, with the same 
Jastrow factor for all states: 



k 



V kt O ki ) 2 )i 



(20) 



where the local energy and the functions O k i are defined 
as before but refer now to the different trial states, 



E 



(i) 



and O k i 



(21) 



and (-)i denotes the average with respect to l^il 2 . Since 
the optimization of orbitals is expressed as optimization 
of CI coefficients in the external space, all parameters 
become state-specific and the minimization of \ 2 leads 
to a different set of linear equations for each state. The 
corrections V k i from these equations are then combined 
with the weights to yield a new single approximate 
effective Hamiltonian. 

After this averaged effective Hamiltonian has been di- 
agonalized, a set of natural orbitals common to all states 
is obtained from the averaged single particle density ma- 
trix 



PSA 



(22) 



ieA 



As a result of this last step, the corrections Vm cor- 
responding to external excitations do not vanish upon 
convergence and, as in the SA-MCSCF, only the aver- 
aged energy will become stationary with respect to or- 
bital variations. 

Since our procedure imposes orthogonality only among 
the determinantal parts the full wave functions 
will in general not be orthogonal due to the presence of 
the Jastrow factor. Therefore, the generalized variational 
theorem providing lower bounds for the excited state en- 
ergies can be violated. 



D. The stochastic reconfiguration approach to 
wave function optimization 

We briefly review the so-called stochastic reconfigu- 
ration (SR) approach22i2£ of Sorella et al. and discuss 
similarities and differences with the EFP method. 

In the SR method, an improved state is obtained by 
applying the operator A — TL to the current trial state 



#(n+i) = (A - U)^ {n) 



(23) 



where A is an energy shift which controls the rate of 
convergence. 

If there is a parameterized ansatz for \J/W with param- 
eters a k ,k = 1 . . .p, a linear variation with respect to all 
parameters can be written as 



(24) 



fe=0 



where the quantities O k with k > are the logarithmic 
derivatives of \&' n ) with respect to the parameter a k 



aa k 



(25) 



and S k denotes a change in the parameter a k . The oper- 
ator Oq is the identity and the associated parameter ao 
corresponds to an overall scaling of the wave function. 

The central idea of the SR approach is to apply A — H 
to the current trial wave function and project the result 
onto the space defined by the parameterization. This 
leads to the conditions 

(* (n) |O fe (A - H)\¥ n) ) = (*(")|O fe |* (n+1) ). (26) 

After inserting Eq.|^]for \£r( n+1 ) and replacing with 
^t, the following equations are obtained for k = . . . m: 



{O k (A-E L ))=J2^{OkO l ). 



(27) 



The k = equation yields the scaling Sq in terms of A: 
5 = A-E-Y,Si(O l ). (28) 



;>o 



After substitution of A, the equations with k > become 
E(O k )-(O k E L ) =J2H(OkOi)-(O k )(O l )), (29) 



which can be rewritten as 

- (AO k AE) =Y / S l (AOiAO k ) 



(30) 



z>o 



These equations coincide, apart from the minus sign, 
with the working equations of the EFP method (Eq. [§J| . 
However, they describe different quantities: changes in 
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the wave function (5i) in the SR method, and changes to 
the Hamiltonian (Vi) in the EFP method. 

To investigate the connection between the EFP and the 
SR method, let us suppose that the new eigenstates of the 
EFP effective Hamiltonian are not obtained through di- 
agonalization but approximately to first order in the per- 
turbation given by the corrections VJ. Then, the change 
of the trial state with respect to the fc-th eigenstate 
of the old effective Hamiltonian is given by 

*T = J$0 J (*0 - E ^! E ®k) > ( 31 ) 

while, in the SR method, the same quantity is obtained 
as 

*t = J®o J (&o + |j*fc) , (32) 

where the parameter 5k corresponds to a variation with 
respect to the eigenstate Therefore, the main differ- 
ence between the two methods seems to be a parameter 
specific scaling in the EFP as opposed to a global one 
in the SR method. We found this difference to affect the 
convergence rate of the methods sometimes considerably, 
as will be shown in Sec. IIVI 

As pointed out by Sorella et al, there is a relation- 
ship between the choice of the energy shift parameter 
A and the amount of sampling to determine the correc- 
tions 5k in each iteration. The more accurately these are 
sampled, the smaller one can make A, achieving faster 
convergence. However, as shown in Appendix [HI a sim- 
ple consideration of the convergence behavior suggests 
that there exists a critical A c below which convergence 
can not be achieved, a finding which is in agreement with 
our numerical results. 



III. COMPUTATIONAL DETAILS 

The vertical excitations of ethene are computed using 
the experimental ground state geometry which is of T)2h 
symmetry (i? C H = 1.0861, R C c = 1.3391 and ZHCH = 
117.6°)2i. The molecule is placed in the yz- plane with 
the molecular axis along the z direction. 

The carbon Is electrons are replaced by a norm- 
conserving s-non-local pseudopotential generated in an 
all-electron Hartree-Fock calculation for the carbon 
atom. The potential of the hydrogen atom is softened 
by removing the Coulombic divergence. 

The Gaussian basis sets are optimized for our soft 
pseudopotentials and augmented with polarization and 
diffuse functions. The calculations are performed with 
two different basis sets. Basis (A) is a contracted 
(12sl2p2d) /[5s5p2d\ basis with the most diffuse exponent 
being 0.02. Basis (B) consists of basis (A) augmented 
with two more diffuse s and p functions with exponents 
0.005 and 0.002. 



The Hartree-Fock and CASSCF calculations are per- 
formed with the program GAMESS(US)^. In all SA- 
MCSCF calculations, equal weights are employed for the 
two states. The program package CHAMR 3 ^. is used for 
the quantum Monte Carlo calculations. Different Jas- 
trow factors are used to describe the correlation with 
a hydrogen and a carbon atom, and their parameters 
are optimized within QMC using the variance minimiza- 
tion method^. We employ both a '2-body' Jastrow fac- 
tor consisting of electron-electron and electron-nucleus 
terms, and a '3-body' Jastrow factor where additionally 
electron-electron-nucleus terms are included^. An imag- 
inary time step of 0.075 H" 1 is used in the DMC calcu- 
lations. 

Singular value decomposition with a threshold of 
0.0001 is used for inverting the matrix (AO/AOfe) in 
Eq. El We always find a large gap in the spectrum of 
this matrix so that the results are not sensitive to the 
precise value of the threshold. 



IV. RESULTS AND DISCUSSION 

All traditional quantum chemistry as well as QMC 
techniques rely on the fact that their reference or trial 
wave function captures the essential nature of the state 
under consideration, whose description is then refined by 
including dynamic correlation. This precisely fails in the 
case of the valence l 1 B lu state of ethene since explicit 
inclusion of dynamic correlation already in the reference 
seems necessary to avoid mixing with Rydberg states. 
Moreover, the problem worsens with increasing single- 
particle basis set, especially with the addition of more 
diffuse functions, since a simple reference wave function 
yields a state of even more Rydberg nature. 

An additional complication with this state is posed by 
the occupied a orbitals which respond to the 7r — > n* exci- 
tation and cannot be treated as frozen. Although ethene 
seems to be a very special case, similar problems are 
expected to be present in many photoactive molecules, 
possibly in a milder but also less clear-cut way. For all 
these reasons, the l 1 Bi ll state of ethene represents a very 
stringent test for methods which aspire to provide highly 
accurate excitation energies, and is therefore chosen here 
as an illustrative example to demonstrate the effective- 
ness of our EFP optimization method. 

Although the ^Bi^ state of ethene is an ideal test 
case for sophisticated correlation methods, a down-side 
must be mentioned: it is very difficult to extract reliable 
estimates of this vertical excitation from gas phase ex- 
periments since ethene starts immediately to twist upon 
photoexcitation. Current interpretations of the measured 
data seem to indicate 7.7 eV as a lower bound to the 
l 1 Bi ll vertical excitation energy of ethen o 35 ! 36 . Theo- 
retical excitation energies obtained with contemporary 
quantum chemistry methods vary between 7.69 eV from 
MRCSi and 8.4 eV from CASPT2i&. 

In Table H we present the VMC and fixed-node DMC 
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TABLE I: Total VMC and DMC energies in Hartree and spatial extent of the wave function {X 2 } in a 2 , for the states l^Ag,, 
l 1 Biu, 2 1 A 9 , 2 1 Bi„ and 3 1 Bi tl of ethene. Different combinations of basis sets and active spaces are used. A 2-body Jastrow 
factor is employed unless indicated as '3body-J'. The total number of occupied orbitals in the reference is listed for each wave 
function type, together with the number of optimized orbitals and variational parameters. The DMC excitation energies in 
eV are computed with respect to the ground state energy obtained with the same basis and Jastrow factor as in the excited 
state. We do not report the ground state DMC energies for basis (B) for reoptimized orbitals or a 3-body Jastrow factor since 
the same behavior is observed as for basis (A). The numbers in parentheses are the statistical errors. See text for a detailed 
explanation. 



state 


basis 


wave function 


occupied orb. 


optimized orb. 


parameters 




(X 2 ) 




AE (eV) 


1 1 A 9 


A 


HF 


6 





- 


-13.6744(5) 


12 


-13.7194(4) 


- 










6 


101 


-13.6797(5) 


12 


-13.7204(4) 


- 






3body-J 




6 


101 


-13.6935(5) 


11 


-13.7203(4) 


- 




B 


HF 


6 





- 


-13.6737(5) 


12 


-13.7191(4) 


- 


l 1 Bi„ 


A 


CAS 2-2 


7 





- 


-13.3502(5) 


37 


-13.4095(4) 


8.45(2) 








•j a 





- 


-13.3321(5) 


12 


-13.4116(4) 


8.40(2) 










2 


17 


-13.3694(5) 


20 


-13.4245(4) 


8.05(2) 










7 


201 


-13.3738(5) 


20 


-13.4257(4) 


8.02(2) 






3body-J 




7 


201 


-13.4011(5) 


19 


-13.4292(4) 


7.92(2) 






CAS 6-6 


9 





- 


-13.3546(5) 


32 


-13.4134(4) 


8.35(2) 










9 


858 


-13.3803(5) 


18 


-13.4280(4) 


7.95(2) 






3body-J 




9 


858 


-13.4008(5) 


19 


-13.4285(5) 


7.94(2) 




B 


CAS 2-2 


7 







-13.3517(5) 


52 


-13.4091(4) 


8.47(2) 










7 


233 


-13.3739(5) 


22 


-13.4245(4) 


8.05(2) 






3body-J 




7 


233 


-13.4004(5) 


24 


-13.4289(4) 


7.93(2) 


2% 


A 


CAS 2-2 


7 







-13.3687(5) 


52 


-13.4279(4) 


7.96(2) 










7 


293 


-13.3708(5) 


46 


-13.4129(4) 


8.36(2) 






3body-J 




7 


293 


-13.3868(4) 


49 


-13.4129(4) 


8.36(2) 


2 1 Bi tl 


A 


CAS 2-3 


8 







-13.3058(5) 


42 


-13.3877(4) 


9.05(2) 










8 


401 


-13.3136(5) 


49 


-13.3699(4) 


9.53(2) 






3body-J 




8 


401 


-13.3439(4) 


52 


-13.3758(4) 


9.37(2) 


^Biu (SA)» 




CAS 2-3 


8 


8 


401 


-13.3718(5) 


22 


-13.4233(4) 


8.08(2) 






3body-J 




8 


401 


-13.3994(5) 


20 


-13.4279(4) 


7.95(2) 


l 3 Bi u 


A 


CAS 2-2 


7 


7 


293 


-13.5118(5) 


12 


-13.5553(4) 


4.49(2) 



a Two-determinant wave function with triplet orbitals. 
b This state is obtained in a state-average (SA) calculation with the corresponding 2 1 Bi tl state. 



TABLE II: Comparison of vertical excitation energies of 
ethene from fixed-node DMC with experiment and other the- 
oretical studies. All excitation energies are in eV. 



State 


DMC 


Exp. 


MRCI 




7.93(2) 


>7.7 a 


7.69 b , 7.8 ... 7.9 C , 7.96 d 


l 3 Biu 


4.49(2) 4.50(3) e 


4.36 d , 


4.49 d 


2 l A g 


8.36(2) 


8.29 d 


8.21 d 


2 -Bin 


9.37(2) 


9.33 d 


9.17 d 



a Ref. GH b Ref. EJ c Ref. d d Ref. Q e Ref. 



energies for the ground state 1 1 A S and the l 1 Bi„ state 
of ethene. We also list the VMC expectation value of the 
spread of the wave functions in the direction perpendicu- 
lar to the molecular plane, (X 2 ) = ( , Ft| Yli x ?I^t}, since 



this quantity has proven useful to differentiate between 
the valence and Rydberg nature of a state. Additionally, 
we show results for the higher singlet states 2 1 A g and 
2 1 Bi„ and the l 3 Bi„ triplet state. 

For the l 1 Bi u state, we investigate how the QMC exci- 
tation energy varies when going from basis set (A) to the 
more diffuse basis set (B), and whether our EFP opti- 
mization method is able in both cases to correct the Ry- 
dberg character of the starting MCSCF reference used in 
the Jastrow-Slater wave function. Since the l 1 Bi u state 
is expected to have a pronounced HOMO-LUMO charac- 
ter, we first employ a simple two-determinant wave func- 
tion, corresponding to a 2 electron in 2 orbitals CASSCF 
wave function (denoted by 'CAS 2-2'). To improve on 
possible deficiencies of this description, we also consider 
a 6 electron in 6 orbitals CASSCF wave function (de- 
noted by 'CAS 6-6') where the active space consists of 
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the orbitals lb2 g , 1&3«, 2a gi \b^ g , 2&2u, and 2bi u . Since 
optimizing both CI coefficients and orbitals for this ac- 
tive space would result in more than 10 5 CSF's which 
is not feasible with our present implementation, we first 
optimize the CI coefficients with respect to a fixed set 
of orbitals. Upon convergence, we apply a threshold of 
0.01 to the coefficients of the CSF's and augment the re- 
sulting truncated CI expansion with single excitations in 
order to relax the orbitals. In this last step, single ex- 
citations from 'core' orbitals, i.e. occupied orbitals that 
were not included in the CAS, can be considered as well. 
The number of virtual orbitals for single excitations is 
dictated by technical limitations: while, for the CAS 2-2 
wave function, all virtual orbitals can be used, in the case 
of the truncated CAS 6-6 wave function, we include the 
lowest 50 virtual orbitals. An illustration of the two-step 
procedure is given in Fig. ^ after the initial optimiza- 
tion of the CI coefficients, the orbitals in the reference 
comprised by the relevant CSF's are relaxed. This two- 
step approach is justified by the very different role of the 
reference wave function in QMC compared to traditional 
quantum chemistry methods: a smaller number of de- 
terminants is needed in a Jastrow-Slater wave function 
since the reference does not define the available excita- 
tion space for the description of dynamical correlation as 
it is the case for a method like MRCI. Finally, we inves- 
tigate the effect on the excitation energy as an increasing 
number of occupied orbitals are reoptimized, and the im- 
pact of using a 3-body versus a 2-body Jastrow factor. 

The results in Table UJ for the l 1 Bi„ state can be sum- 
marized as follows. The initial MCSCF Jastrow-Slater 
wave functions have a substantial Rydberg contribution 
which, as expected, increases when the basis set becomes 
more diffuse. The spread (A 2 ) in the direction perpen- 
dicular to the molecular plane increases from 37 a 2 , with 
basis (A) to 52 a 2 , with basis (B) for the CAS 2-2 wave 
function. The resulting fixed-node DMC excitation ener- 
gies are about 8.5 eV, thus much higher than what was 
found in other benchmark ab initio calculations. They 
roughly agree with CASPT2 results which are believed 
to suffer from the same deficiency. 

Optimizing the orbitals in the presence of the Jastrow 
factor, i.e. including the feedback of dynamic correlation 
on the reference, reduces the spread of the wave function 
to 18-22 do, slightly larger than the ground state value. 
The DMC excitation energies are lowered by as much as 
0.5-0.6 eV. After optimization, the results obtained with 
basis A and basis B are indistinguishable. 

For CAS 2-2 wave functions, the most substantial re- 
duction in the VMC and DMC energies and in (A 2 ) is 
obtained when the active orbitals are optimized in the 
presence of the Jastrow factor. If a orbitals are included 
in the optimization, further smaller improvements are 
gained at the variational level while the DMC energies 
remain rather insensitive. Using a CAS 6-6 wave func- 
tion instead of the two-determinant wave function yields 
lower VMC energies but, for fixed-node DMC, the effect 
of the multi-configuration trial wave function is not very 



large. In fact, after the optimization, there remain only 
10 determinants with CI coefficients above our thresh- 
old of 0.01 while, in the initial CAS 6-6 wave function 
from MCSCF, 36 determinants meet the same criterion. 
When a 3-body Jastrow factor is employed, the differ- 
ence between a CAS 6-6 and a CAS 2-2 energy becomes 
negligible also at the VMC level. In general, the use 
of a 3-body instead of a 2-body Jastrow factor always 
gives lower VMC energies and lower excited-state DMC 
energies^ while the DMC energies of the ground state 
are rather unaffected. The resulting improvement on the 
DMC excitation energy is visible only when a CAS 2-2 
Jastrow-Slater wave function is used to localize the pseu- 
dopotentiali. 



-13.36, 

-13.37 

H -13.38 

5,-13.39 
m 

-13.40 
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FIG. 1: Example of the convergence of the VMC energy and 
the spatial extent (X 2 ) of the wave function in the EFP opti- 
mization of a CAS 6-6 wave function with basis (A). Initially, 
only the CI coefficients are being optimized. After the first 10 
iterations, a threshold of 0.01 is applied to the coefficients of 
the CSF's and only single excitations on top of this truncated 
CI wave functions are considered to optimize all the occupied 
orbitals. The statistical error is smaller than the size of the 
symbols. See the text for more details. 

The reduction of the spatial extent of the trial wave 
function after reoptimization to values around 20 arj 2 is 
in qualitative agreement with the findings of Buenker et 
al^m and Lischka et alw^ although we observe a consider- 
ably larger spread of the values depending on the details 
of the optimization. As Fig.njillustrates, the convergence 
of (A 2 ) with the number of steps in the EFP optimiza- 
tion is much slower than the convergence of the energy 
and harder to establish. Interestingly, the EFP optimiza- 
tion appears to be sensitive to quantities other than the 
energy: even though the energy is practically converged, 
additional optimization steps yield a systematic lowering 
of (A 2 ). 

Our final DMC results are summarized and compared 
with experiment and other theoretical studies in Table ITU 
The best estimate for the excitation energy of the l 1 Bi u 
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state is 7.92 ± 0.02 eV, in good agreement with most 
other ab initio calculations^* 1 *!. However, it is higher by 
about 0.23 eV than the value of 7.69 eV recently obtained 
by Lischka et al. using a sophisticated MRCI approach 
with up to 12 electrons in 12 orbitals reference spaces. It 
should be noted that such MRCI calculations suffer from 
a size-consistency problem which can be approximately 
corrected using a term proposed by Davidson (see e.g. 
Ref . l2l|) . For the present example, the Davidson correc- 
tion amounts to about 0.1 eV and Bucnker and Krebs 
questioned its trustworthiness^ 9 -. 

The MRCI results of Buenker and Krebs were partly 
obtained using MCSCF triplet orbitals which are as com- 
pact as the ground state orbitals. Following this estab- 
lished strategy, recent DMC calculations also employed 
triplet orbitals in trial wave functions for singlet states^. 
However, when we use this recipe for ethene, the corre- 
sponding fixed-node DMC value is 8.4±0.02 eV, which is 
comparable with the value obtained from singlet MCSCF 
orbitals. 

In the case of the 2 1 A S and 2 1 Bi u states, we employ 
our state-averaged EFP optimization approach because 
these states are not the lowest in their irreducible repre- 
sentation and, therefore, have to be optimized in the pres- 
ence of the 1 1 A S and l 1 !?^ states, respectively. In DMC, 
we rely instead on the fixed-node approximation to pre- 
vent a collapse to the lower states (see AppendixEJl ■ The 
EFP optimization of multiple states is found to be quite 
stable: both VMC and DMC energies of the l 1 i?i u state 
taken from the state-averaged optimization (denoted by 
'SA' in Table D) are in agreement with the correspond- 
ing values obtained by optimizing the lower state alone. 
Moreover, upon optimization, the DMC excited state en- 
ergies are substantially higher than those obtained us- 
ing MCSCF orbitals which give excitation energies lower 
than the experimental values by 0.33 eV for the 2 : A 5 
and 0.48 eV for the 2 1 B\ U state. This can possibly be 
understood from the fact that valence- Rydberg mixing in 
the initial MCSCF wave function raises the lower state 
and lowers the higher one. As in the case of the l 1 i3i u 
state, the best excitation energy for the 2 1 B\ U state is ob- 
tained when a 3-body Jastrow factor is employed, while 
the 2 1 A g excitation is insensitive to the change from a 2- 
body to a 3-body Jastrow factor. Our final DMC excita- 
tion energies of 8.36±0.02 eV for 2 1 A g and 9.37±0.02 eV 
for 2 l B\ u are only slightly higher than the correspond- 
ing experimental values^ of 8.29 eV and 9.33 eV, respec- 
tively. 

For the l 3 Bi u state we obtain an excitation energy of 
4.49±0.02 eV which agrees with the recent DMC resuh^l 
of El Akraminc et al., the experimental value quoted in 
Refii4 and the MRCI calculation of Krebs and Bucnker 19 . 
The higher MRCI value of 4.61eV reported by Gemein 
and Pcycrimhoff 18 seems to result from using a somewhat 
different geometry. 

Finally, we examine the differences between our EFP 
approach and the stochastic reconfiguration method by 
Sorella et al.. As explained in Sec. Ill Dl the SR method 




VMC Energy 



o SR A = - 10 
□ SR A = -6 
o EFP 
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— // — 

Iteration 
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FIG. 2: Convergence of the VMC energy and the spatial ex- 
tent (X 2 ) of the ethene I 1 Bi u state with the iteration number 
in the optimization. The EFP and the SR method with two 
different values of the parameter A are used to optimize the 
two active orbitals in a two-determinant wave function with 
basis (A). For values of A lower than -10 a.u., the SR opti- 
mization is unstable. The statistical error is smaller than the 
size of the symbols. See text for more details. 



applied to the optimization of the determinantal part of 
the trial wave function can be considered as a simpler 
variant of the EFP approach that should lead to the same 
solution upon convergence but may have different conver- 
gence properties. In Fig. the different performance of 
the two approaches is illustrated with the optimization 
of the two active orbitals in a two-determinant basis- (A) 
wave function of the l 1 Bi 11 state of ethene. Both calcu- 
lations are performed using the same amount of Monte 
Carlo sampling per iteration. While the EFP method is 
parameter-free, the convergence of the SR approach de- 
pends on the choice of the parameter A. We find the 
optimal value of A leading to the fastest convergence to 
be about —10 a.u., a value which approximately agrees 
with our simple estimate of the critical value A c given 
in Appendix El Values of A higher than A c yield slower 
convergence as can be seen in Fig. EJ On the other hand, 
decreasing A below A c results in the divergence of the en- 
ergy since the contributions of some excited states in the 
wave function start to be amplified instead of damped. 
While the EFP optimization is essentially converged after 
about 10 iterations, the SR approach even at the optimal 
value of A = —10 a.u. takes more than 100 iterations to 
converge. As explained in Appendix El the critical value 
A c is related to the energy spread inside the variational 
space which, in the case of orbital variations, is the space 
spanned by the single excitations from the reference con- 
figurations into the set of admitted virtual orbitals. For 
basis (A) and a two-determinant reference, if we use all 
virtual orbitals of proper symmetry, this spread amounts 
to about 6 a.u. If we reduce the energy spread by exclud- 
ing high-lying virtual orbitals which are not expected to 
contribute significantly to the optimization, A can be fur- 
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ther decreased, yielding a faster convergence. However, 
the convergence of the SR approach is expected to be- 
come slower as the size of the system increases since the 
energy spread of the variational space will increase as 
well. 



V. CONCLUSIONS 

We extended the energy fluctuation potential (EFP) 
method to simultaneously optimize the orbitals and the 
CI coefficients in Jastrow-Slater wave functions via a 
super-CI approach, and to treat state-averaging for the 
optimization of multiple states of the same symmetry. 
With these additions, the method becomes a useful and 
effective tool to optimize trial wave functions for both 
ground and excited states, which can then be used in 
fixed-node diffusion Monte Carlo. 

As illustrative examples, we considered several vertical 
excitations of ethene, in particular the difficult valence 
^Biu state. For this state, the EFP approach leads to 
strongly improved trial wave functions upon the starting 
Rydberg-like MCSCF reference, as can especially be seen 
from the reduced spatial spread of the optimized wave 
function. All results seem in reasonable agreement with 
experiment and other quantum chemical methods. In 
contrast to techniques like MRCI and CASPT2, it is not 
necessary to use multi-configuration trial wave functions 
once the orbitals are optimized in the presence of the 
Jastrow factor. Also, using triplet orbitals for singlet 
states does not appear to be a reliable recipe for fixed- 
node DMC. 

Our EFP optimization of multiple states was employed 
for the two Rydberg-like states 2 1 A g and 2 1 B lu of ethene, 
which are not the lowest ones in their irreducible rep- 
resentations. The optimization procedure is stable and 
yields substantially higher DMC excited state energies 
than the ones obtained with MCSCF orbitals, bringing 
the corresponding excitation energies in better agreement 
with the experimental values. 

In addition, we also compared the EFP method to 
the stochastic reconfiguration method of Sorella et al~. 
These two methods stem from quite different theoreti- 
cal backgrounds but, for practical purposes, the EFP ap- 
proach can be regarded as a more sophisticated variant of 
the SR scheme. In general, the EFP method shows much 
faster convergence than the SR approach. On the other 
hand, the SR technique can be applied to the optimiza- 
tion of arbitrarily parameterized wave functions while the 
EFP method is always based on a suitable factorization 
of the wave function. 
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APPENDIX A: ORBITAL TRANSFORMATION 

Let U be the orthogonal transformation between the 
current orbitals <f> and the corresponding natural orbitals 
4> nat ref of the reference wave function 

M 

= £t^. (Al) 

i=i 

After updating and rediagonalizing the super-CI Hamil- 
tonian (i.e. reference plus single excitations), the natural 
orbitals of this wave function are given by the transfor- 
mation V: 

N 

^natext = £ V . j(j) . (A2 ) 

i=i 

Here, M and N (N 3> M ) denote the dimensions of the 
reference space and the super-CI space, respectively. The 
new reference orbitals <fi are now obtained as 

M N 

■'■ Y^l 1 - V -^ *=1-..M, (A3) 

j=i k=i 

where the orthogonality of U (U^ 1 = U T ) was used. 
Upon convergence, this ensures that the natural orbitals 
of the reference wave function coincide with the corre- 
sponding subset of the natural orbitals of the super-CI 
expansion. 

As the new virtual orbitals [4>i,i = M + 1 . . . N) one 
can use the orbitals <f>f at cxt , i — M + 1 . . . N since they 
are orthogonal to the occupied ones. Alternatively, one 
can obtain a new set of virtual orbitals by explicit or- 
thogonalization of the orbitals of the previous iteration 
to the new occupied orbitals. 

It should be noted that in the case of a complete ref- 
erence space (CAS) the natural orbitals of the super-CI 
expansion could be used instead of the orbitals from the 
transformation above, since the CAS wave function is 
invariant with respect to transformations among the ac- 
tive orbitals. However, subsequent truncation of the wave 
function (omission of determinants with coefficients be- 
low some threshold) could lead to slightly different results 
depending on the orbital transformation employed. 

APPENDIX B: CRITICAL SCALING 
PARAMETER OF THE SR METHOD 

Let us suppose that the M eigenstates of some Hamil- 
tonian Ti are known 

H*i=Ei9 t , i = 0,...,M. (Bl) 
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An approximation \l/ to the state ^q can be written as 

i>0 

The variations with respect to the states i > are intro- 
duced as 

*' = * + £ = * ° + £ • ( B3 ) 

i>0 i>0 

Acting with the operator A — Ti on 'J yields 
* = (A-W)* 

= (A - £ )*o + C ^ A ~ ■ ( B4 ) 

Therefore, after dividing by A — Eq and equating ^' with 
'J, we obtain that the corrections on are given by 



on 



A - .En 



(B5) 



Here, there is no need for a projection since applying 
A — H does not give rise to contributions outside the 
variational space. The method converges if 



M<1 Vi, 



(B6) 



where on = would give convergence in one step while 
ratios closer to ±1 yield slower convergence. This trans- 
lates into the condition A > A c = {Eq + Em)/2- At 
the critical value A c , the contribution due to 'J'i will be 
suppressed at the rate 



Em — Eq + 2(Eq — Ej) 
Em — Eq 



(B7) 



which becomes smaller as the energy spread in the vari- 
ational subspace Em — Eq increases. 



APPENDIX C: FIXED NODE DIFFUSION 
MONTE CARLO AND EXCITED STATES 

In diffusion quantum Monte Carlo (DMC), the imag- 
inary time evolution operator is used to stochastically 
project a trial wave function onto the lowest energy 
state $o of the system to which is not orthogonal: 



* = lim exp {-tH} *t 



(CI) 



For fermions, the so-called fixed-node approximation is 
generally introduced in order to prevent collapse to the 
bosonic solution: the result of the projection is con- 
strained to have the same nodal surface as a given trial 
wave function ^t- This is equivalent to separately solv- 
ing the Schrodinger equation in the regions of constant 
sign (nodal pockets) of \&t subject to Dirichlet bound- 
ary conditions. The nodal pockets of a ground state wave 
function can be shown to be equivalent, and good ground 
state trial wave functions seem to approximate this prop- 
erty well. 

In the same way as a trial wave function can be used 
within the fixed-node approximation to prevent collapse 
to the energetically lower bosonic state, it can also be em- 
ployed for an excited state calculation to prevent collapse 
to lower states. If the excited state trial wave function 
has the exact nodes, the exact excited state energy is re- 
covered. However, with approximate nodes, the method 
is not variational except for the lowest state of each one- 
dimensional irreducible representation of the pointgroup 
of the molecule 38 . 
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